function bits = HelperGPSNAVDataEncode(cfg)
%HelperGPSNAVDataEncode GPS CNAV/LNAV Navigation data bits generation
%
%   Note: This is a helper and its API and/or functionality may change
%   in subsequent releases.
%
%   BITS = HelperGPSNAVDataEncode(CFG) generates encoded bit stream BITS
%   which contains data bits for GPS legacy navigation (LNAV) or civil
%   Navigation (CNAV) data transmission. CFG is a configuration object of
%   type <a href="matlab:help('HelperGPSNavigationConfig')">HelperGPSNavigationConfig</a>. The data bits are generated from the
%   fields that go into legacy navigation (LNAV) GPS data. Each of these
%   fields are available in the config object, CFG. The frame structure of
%   LNAV data, CNAV data and the properties that are encoded into the frame
%   are given in Appendix II and Appendix III of IS-GPS-200L [1]
%   respectively.
%
%   References:
%    [1] IS-GPS-200L. "NAVSTAR GPS Space Segment/Navigation User Segment
%        Interfaces." GPS Enterprise Space & Missile Systems Center (SMC) -
%        LAAFB, May 14, 2020.
%
%   See also HelperGPSNavigationConfig, HelperGPSCEIConfig,
%   HelperGPSBasebandWaveform.

%   Copyright 2021-2022 The MathWorks, Inc.

if strcmp(cfg.DataType,'CNAV')
    bits = GPSCNAVDataEncode(cfg);
else % strcmp(cfg.DataType,'LNAV')
    bits = GPSLNAVDataEncode(cfg);
end
end

function navData = GPSLNAVDataEncode(cfg)
%GPSLNAVDataEncode Navigation data bits generation for LNAV GPS
%
%   LNAVBITS = GPSLNAVDataEncode(CFG) generates encoded bit stream
%   LNAVBITS which contains data bits for legacy GPS transmission. CFG is
%   a configuration object of type <a href="matlab:help('HelperGPSNavigationConfig')">HelperGPSNavigationConfig</a>. The data bits
%   are generated from the fields that go into legacy navigation (LNAV) GPS
%   data. Each of these fields are available in the config object, CFG.
%   The frame structure of LNAV data and the properties that are encoded
%   into the frame are given in Appendix II of IS-GPS-200L [1].

subframeLength = 300; % Bits
NumSubframes   = 5;
frameIndices = cfg.FrameIndices;
numFrames = length(frameIndices);
navData = zeros(numFrames*NumSubframes*subframeLength, 1);

for iFrame = 1:numFrames
    cfg.PageID = modlocal(cfg.FrameIndices(iFrame), 25);
    for iSubframe = 1:NumSubframes
        cfg.SubframeID = iSubframe; % This information is encoded in every subframe. Hence, need to be updated for every subframe
        switch iSubframe
            case 1
                subframe = subframe1(cfg);
            case 2
                subframe = subframe2(cfg);
            case 3
                subframe = subframe3(cfg);
            case 4
                subframe = subframe4(cfg);
            case 5
                subframe = subframe5(cfg);
        end
        cfg.HOWTOW = cfg.HOWTOW + 1; % The TOW in HOW increments by 1 for every 6 sec.
        startIDX = (iFrame-1)*(1500) + 1 + (iSubframe-1)*300;
        endIDX = (iFrame-1)*(1500) + iSubframe*300;
        navData(startIDX:endIDX) = subframe;
    end
end
end

function word = HandoverWord(howdata)
%HANDOVERWORD generates the handover word (HOW)
%   WORD = HandoverWord(HOWDATA) generates the 30 bits for the handover
%   word according to the parameters given in the structure HOWDATA.
%   Following are the parameters in the structure HOWDATA:
%
%   HOWTOW           - time of week in the handover word. LSB is 6 sec
%   CEIDataSet.AlertFlag     - Alert flag
%   AntiSpoofFlag    - Anti-spoofing flag. If set, then anti-spoofing is on
%   SubframeID       - Subframe ID (1 to 5)

towBin = HelperGPSNum2BitsVector(howdata.HOWTOW, 17, 1);
bit18Alert = howdata.CEIDataSet.AlertFlag;
bit19AS = howdata.AntiSpoofFlag;
subframeIDBin = HelperGPSNum2BitsVector(howdata.SubframeID, 3, 1);
data = [towBin; bit18Alert; bit19AS; subframeIDBin];
wordNumber = 2;
word = gpsLNavWordEnc(data,wordNumber);
end

function reservedWord = ReservedWord(wordNumber)
% Generate a word whose 24 bits data elements are all zeros and have valid
% parity bits. As parity calculation depends on previous word, parity is not
% always zeros even if 24 bits in the word is zero.

data = zeros(24, 1);
reservedWord = gpsLNavWordEnc(data,wordNumber);
end

function subframe = subframe1(opts)
% Generate 300 bits of GPS LNAV subframe1 data from the parameters given in
% 'opts' input structure

word1 = TLMWord(opts); % Having the default values of TLM word for speed
word2 = HandoverWord(opts);

wordNumber = 3;

weekNumber = mod(opts.CEIDataSet.WeekNumber,1024);
weekNumberBin = HelperGPSNum2BitsVector(weekNumber,10, 1);

if strcmp(opts.CodesOnL2,'P-code')
    codesOnL2Bin = [0;1];
else % C/A-code
    codesOnL2Bin = [1;0];
end

uraVal = opts.CEIDataSet.URAID;
uraValBin = HelperGPSNum2BitsVector(uraVal, 4, 1);

svHealthBin = HelperGPSNum2BitsVector(opts.CEIDataSet.SVHealth, 6, 1);

iodc = opts.CEIDataSet.IssueOfDataClock;
iodcBin = HelperGPSNum2BitsVector(iodc, 10, 1);

data = [weekNumberBin; codesOnL2Bin; uraValBin; svHealthBin; iodcBin(1:2)];
word3 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 4;
data = [opts.L2PDataFlag; zeros(23, 1)];
word4 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 5;
word5 = ReservedWord(wordNumber);

wordNumber = 6;
word6 = ReservedWord(wordNumber);

wordNumber = 7;

T_GDBin = HelperGPSNum2BitsVector(opts.CEIDataSet.GroupDelayDifferential, 8, 2^-31);

data = [zeros(16, 1); T_GDBin];
word7 = gpsLNavWordEnc(data, wordNumber);

wordNumber = 8;

tOC = opts.CEIDataSet.ReferenceTimeOfClock;
tOCBin = HelperGPSNum2BitsVector(tOC, 16, 2^4);

data = [iodcBin(3:10); tOCBin];
word8 = gpsLNavWordEnc(data,wordNumber);

af2 = opts.CEIDataSet.SVClockCorrectionCoefficients(3);
af1 = opts.CEIDataSet.SVClockCorrectionCoefficients(2);
af0 = opts.CEIDataSet.SVClockCorrectionCoefficients(1);

wordNumber = 9;

scaleFactor = 2^(-55);
af2Bin = HelperGPSNum2BitsVector(af2, 8, scaleFactor);

scaleFactor = 2^(-43);
af1NumBits = 16;
af1Bin = HelperGPSNum2BitsVector(af1, af1NumBits, scaleFactor);

data = [af2Bin; af1Bin];
word9 = gpsLNavWordEnc(data,wordNumber);

wordNumber=10;

scaleFactor = 2^(-31);
af0NumBits = 22;
af0Bin = HelperGPSNum2BitsVector(af0, af0NumBits, scaleFactor);

data = af0Bin;
word10 = gpsLNavWordEnc(data,wordNumber);

subframe = [word1; word2; word3; word4; word5; word6; word7; word8; word9; word10];

end

function subframe = subframe2(opts)
% Generate 300 bits of GPS LNAV subframe2 data from the parameters given in
% 'opts' and 'ephemeris' input structures

word1 = TLMWord(opts); % Having the default values of TLM word for speed
word2 = HandoverWord(opts);

wordNumber = 3;

iode = opts.CEIDataSet.IssueOfDataEphemeris;
iodeBin = HelperGPSNum2BitsVector(iode,8,1);


crs = opts.CEIDataSet.HarmonicCorrectionTerms(3);
cus = opts.CEIDataSet.HarmonicCorrectionTerms(5);
cuc = opts.CEIDataSet.HarmonicCorrectionTerms(6);
scaleFactor = 2^(-5);
sinHrmcCorrection = HelperGPSNum2BitsVector(crs,16,scaleFactor);

data = [iodeBin; sinHrmcCorrection];
word3 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 4;

meanMotionDiff = opts.CEIDataSet.MeanMotionDifference;
meanMotionDiffBin = HelperGPSNum2BitsVector(meanMotionDiff, 16, 2^(-43));

meanAnom = opts.CEIDataSet.MeanAnomaly;
meanAnomBin = HelperGPSNum2BitsVector(meanAnom, 32, 2^(-31));

data = [meanMotionDiffBin; meanAnomBin(1:8)];
word4 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 5;
data = meanAnomBin(9:32);
word5 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 6;

scaleFactor = 2^(-29);
cucNumBits = 16;
cucBin = HelperGPSNum2BitsVector(cuc, cucNumBits, scaleFactor);

eccen = opts.CEIDataSet.Eccentricity;
scaleFactor = 2^(-33);
eccenNumBits = 32;
eccenBin = HelperGPSNum2BitsVector(eccen, eccenNumBits, scaleFactor);

data = [cucBin; eccenBin(1:8)];
word6 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 7;
data = eccenBin(9:32);
word7 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 8;

scaleFactor = 2^(-29);
cusNumBits = 16;
cusBin = HelperGPSNum2BitsVector(cus, cusNumBits, scaleFactor);

sqRootA = sqrt(opts.CEIDataSet.SemiMajorAxisLength);
scaleFactor = 2^(-19);
sqRootANumBits = 32;
sqRootABin = HelperGPSNum2BitsVector(sqRootA, sqRootANumBits, scaleFactor);

data = [cusBin; sqRootABin(1:8)];
word8 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 9;
data = sqRootABin(9:32);
word9 = gpsLNavWordEnc(data,wordNumber);

wordNumber=10;

toe = opts.CEIDataSet.ReferenceTimeOfEphemeris;
scaleFactor = 2^(4);
toeNumBits = 16;
toeBin = HelperGPSNum2BitsVector(toe, toeNumBits, scaleFactor);

AODO = opts.AgeOfDataOffset;
AODOBin = HelperGPSNum2BitsVector(AODO, 5, 900);

data =[toeBin; opts.CEIDataSet.FitIntervalFlag; AODOBin];
word10 = gpsLNavWordEnc(data,wordNumber);

subframe = [word1; word2; word3; word4; word5; word6; word7; word8; word9; word10];
end


function subframe = subframe3(opts)
% Generate 300 bits of GPS LNAV subframe3 data from the parameters given in
% 'opts' and 'ephemeris' input structures

word1 = TLMWord(opts); % Having the default values of TLM word for speed
word2 = HandoverWord(opts);

wordNumber = 3;

cis = opts.CEIDataSet.HarmonicCorrectionTerms(1);
cic = opts.CEIDataSet.HarmonicCorrectionTerms(2);
crc = opts.CEIDataSet.HarmonicCorrectionTerms(4);
scaleFactor = 2^(-29);
cicNumBits = 16;
cicBin = HelperGPSNum2BitsVector(cic, cicNumBits, scaleFactor);

longOrbitPlane = opts.CEIDataSet.LongitudeOfAscendingNode;
scaleFactor = 2^(-31);
longOrbitPlaneNumBits =32;
longOrbitPlaneBin = HelperGPSNum2BitsVector(longOrbitPlane,longOrbitPlaneNumBits,scaleFactor);

data = [cicBin; longOrbitPlaneBin(1:8)];
word3 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 4;
data = longOrbitPlaneBin(9:32);
word4 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 5;

scaleFactor = 2^(-29);
cisNumBits = 16;
cisBin = HelperGPSNum2BitsVector(cis, cisNumBits, scaleFactor);

refIncAngle = opts.CEIDataSet.Inclination;
scaleFactor = 2^(-31);
refIncAngleNumBits = 32;
refIncAngleBin = HelperGPSNum2BitsVector(refIncAngle, refIncAngleNumBits, scaleFactor);

data = [cisBin; refIncAngleBin(1:8)];
word5 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 6;
data = refIncAngleBin(9:32);
word6 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 7;

scaleFactor = 2^(-5);
crcNumBits = 16;
crcBin = HelperGPSNum2BitsVector(crc,crcNumBits,scaleFactor);

argPerig = opts.CEIDataSet.ArgumentOfPerigee;
argPerigNumBits = 32;
scaleFactor = 2^(-31);
argPerigBin = HelperGPSNum2BitsVector(argPerig,argPerigNumBits,scaleFactor);

data = [crcBin; argPerigBin(1:8)];
word7 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 8;
data = argPerigBin(9:32);
word8 = gpsLNavWordEnc(data,wordNumber);

wordNumber = 9;

rorAscension = opts.CEIDataSet.RateOfRightAscension;
rorAscensionNumBits = 24;
scaleFactor = 2^(-43);
rorAscensionBin = HelperGPSNum2BitsVector(rorAscension,rorAscensionNumBits,scaleFactor);

data = rorAscensionBin;
word9 = gpsLNavWordEnc(data,wordNumber);

wordNumber=10;

iode = opts.CEIDataSet.IssueOfDataEphemeris;
iodeNumBits = 8;
iodeBin = HelperGPSNum2BitsVector(iode, iodeNumBits, 1);

idot = opts.CEIDataSet.InclinationRate;
idotNumBits = 14;
scaleFactor = 2^(-43);
idotBin = HelperGPSNum2BitsVector(idot, idotNumBits, scaleFactor);

data =[iodeBin; idotBin];
word10 = gpsLNavWordEnc(data,wordNumber);

subframe = [word1; word2; word3; word4; word5; word6; word7; word8; word9; word10];
end

function subframe = subframe4(opts)
% Generate 300 bits of GPS LNAV subframe4 data from the parameters given in
% 'opts' and 'almanac' input structures

almanac = opts.Almanac;
word1 = TLMWord(opts); % Having the default values of TLM word for speed
word2 = HandoverWord(opts);

dataid = [0; 1];
subframe4SVID = [57; 25; 26; 27; 28; 57; 29; 30; 31; 32; 57; 62; 52; 53;...
    54; 57; 55; 56; 58; 59; 57; 60; 61; 62; 63];

switch opts.PageID
    case {2, 3, 4, 5, 7, 8, 9, 10} % Almanac data for SV 25 through 32
        prnidx = subframe4SVID(opts.PageID);
        words = almstruct2bits(almanac, prnidx);
    case {1, 6, 11, 12, 14, 15, 16, 19, 20, 21, 22, 23, 24} % Reserved messages
        svid = subframe4SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);

        wordNum = 3;

        data = [dataid; svidBin; zeros(16, 1)];
        word3 = gpsLNavWordEnc(data, wordNum);
        words = [word3; zeros(210, 1)];

        for iword = 4:10
            if iword~=10
                data = zeros(24, 1);
            else
                data = zeros(22, 1);
            end
            words((iword-3)*30+1:(iword-2)*30) = gpsLNavWordEnc(data, iword);
        end
    case 13 % NMCT data
        svid = subframe4SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);

        wordNum = 3;

        ai = opts.NMCTAvailabilityIndicator;
        aiBin = HelperGPSNum2BitsVector(ai, 2, 1);

        nmctBin = NMCT2Bin(opts.NMCTERD);

        data = [dataid; svidBin; aiBin; nmctBin(1, :)'; nmctBin(2,:)'; nmctBin(3, 1:2)'];
        word3 = gpsLNavWordEnc(data, wordNum);

        wordNum = 4;
        data = [nmctBin(3, 3:6)'; nmctBin(4,:)'; nmctBin(5,:)'; nmctBin(6,:)'; nmctBin(7,1:2)'];
        word4 = gpsLNavWordEnc(data, wordNum);

        wordNum = 5;
        data = [nmctBin(7,3:6)'; nmctBin(8,:)'; nmctBin(9,:)'; nmctBin(10,:)'; nmctBin(11,1:2)'];
        word5 = gpsLNavWordEnc(data, wordNum);

        wordNum = 6;
        data = [nmctBin(11,3:6)'; nmctBin(12,:)'; nmctBin(13,:)'; nmctBin(14,:)'; nmctBin(15,1:2)'];
        word6 = gpsLNavWordEnc(data, wordNum);

        wordNum = 7;
        data = [nmctBin(15,3:6)'; nmctBin(16,:)'; nmctBin(17,:)'; nmctBin(18,:)'; nmctBin(19,1:2)'];
        word7 = gpsLNavWordEnc(data, wordNum);

        wordNum = 8;
        data = [nmctBin(19,3:6)'; nmctBin(20,:)'; nmctBin(21,:)'; nmctBin(22,:)'; nmctBin(23,1:2)'];
        word8 = gpsLNavWordEnc(data, wordNum);

        wordNum = 9;
        data = [nmctBin(23,3:6)'; nmctBin(24,:)'; nmctBin(25,:)'; nmctBin(26,:)'; nmctBin(27,1:2)'];
        word9 = gpsLNavWordEnc(data, wordNum);

        wordNum = 10;
        data = [nmctBin(27,3:6)'; nmctBin(28,:)'; nmctBin(29,:)'; nmctBin(30,:)'];
        word10 = gpsLNavWordEnc(data, wordNum);

        words = [word3; word4; word5; word6; word7; word8; word9; word10];
    case 17 % Text message
        svid = subframe4SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);
        numCharInMessage = 22;
        fullText = char(opts.TextMessage(:).');
        if length(fullText) >= numCharInMessage
            text = fullText(1:numCharInMessage);
        else
            text = [fullText, repmat(' ',1, numCharInMessage-length(fullText))];
        end
        TempBits = comm.internal.utilities.de2biBase2LeftMSB((double(text(:))),8)';
        % bits = [bits;GPSText2Bits(text)];

        wordNum = 3;

        data = [dataid; svidBin; reshape(TempBits(:,1:2),[],1)];
        word3 = gpsLNavWordEnc(data, wordNum);

        wordNum = 4;
        data = reshape(TempBits(:,3:5),[],1);
        word4 = gpsLNavWordEnc(data, wordNum);

        wordNum = 5;
        data = reshape(TempBits(:,6:8),[],1);
        word5 = gpsLNavWordEnc(data, wordNum);

        wordNum = 6;
        data = reshape(TempBits(:,9:11),[],1);
        word6 = gpsLNavWordEnc(data, wordNum);

        wordNum = 7;
        data = reshape(TempBits(:,12:14),[],1);
        word7 = gpsLNavWordEnc(data, wordNum);

        wordNum = 8;
        data = reshape(TempBits(:,15:17),[],1);
        word8 = gpsLNavWordEnc(data, wordNum);

        wordNum = 9;
        data = reshape(TempBits(:,18:20),[],1);
        word9 = gpsLNavWordEnc(data, wordNum);

        wordNum = 10;
        data = [reshape(TempBits(:,21:22),[],1); zeros(6,1)];
        word10 = gpsLNavWordEnc(data, wordNum);

        words = [word3; word4; word5; word6; word7; word8; word9; word10];
    case 18 % Ionospheric and UTC data
        ionosphere = opts.Ionosphere;
        utc = opts.UTC;

        wordNum = 3;

        svid = subframe4SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);

        alpha0 = ionosphere.Alpha(1);
        alpha0Bin = HelperGPSNum2BitsVector(alpha0, 8, 2^-30);
        alpha1 = ionosphere.Alpha(2);
        alpha1Bin = HelperGPSNum2BitsVector(alpha1, 8, 2^-27);

        data = [dataid; svidBin; alpha0Bin; alpha1Bin];
        word3 = gpsLNavWordEnc(data, wordNum);

        wordNum = 4;

        alpha2 = ionosphere.Alpha(3);
        alpha2Bin = HelperGPSNum2BitsVector(alpha2, 8, 2^-24);

        alpha3 = ionosphere.Alpha(4);
        alpha3Bin = HelperGPSNum2BitsVector(alpha3, 8, 2^-24);

        beta0 = ionosphere.Beta(1);
        beta0Bin = HelperGPSNum2BitsVector(beta0, 8, 2^11);

        data = [alpha2Bin; alpha3Bin; beta0Bin];
        word4 = gpsLNavWordEnc(data, wordNum);

        wordNum = 5;

        beta1 = ionosphere.Beta(2);
        beta1Bin = HelperGPSNum2BitsVector(beta1, 8, 2^14);

        beta2 = ionosphere.Beta(3);
        beta2Bin = HelperGPSNum2BitsVector(beta2, 8, 2^16);

        beta3 = ionosphere.Beta(4);
        beta3Bin = HelperGPSNum2BitsVector(beta3, 8, 2^16);

        data = [beta1Bin; beta2Bin; beta3Bin];
        word5 = gpsLNavWordEnc(data, wordNum);

        wordNum = 6;

        A1 = utc.UTCTimeCoefficients(2);
        A1Bin = HelperGPSNum2BitsVector(A1, 24, 2^-50);

        word6 = gpsLNavWordEnc(A1Bin, wordNum);

        wordNum = 7;

        A0 = utc.UTCTimeCoefficients(1);
        A0Bin = HelperGPSNum2BitsVector(A0, 32, 2^-30);

        word7 = gpsLNavWordEnc(A0Bin(1:24), wordNum);

        wordNum = 8;

        tot = utc.ReferenceTimeUTCData;
        totBin = HelperGPSNum2BitsVector(tot, 8, 2^12);

        WNt  = utc.TimeDataReferenceWeekNumber;
        WNtBin = HelperGPSNum2BitsVector(WNt, 8, 1);

        data = [A0Bin(25:32); totBin; WNtBin];
        word8 = gpsLNavWordEnc(data, wordNum);

        wordNum = 9;

        DeltLS = utc.PastLeapSecondCount;
        DeltLSBin = HelperGPSNum2BitsVector(DeltLS, 8, 1);

        WNLSF = utc.LeapSecondReferenceWeekNumber;
        WNLSFBin = HelperGPSNum2BitsVector(WNLSF, 8, 1);

        DN = utc.LeapSecondReferenceDayNumber;
        DNBin = HelperGPSNum2BitsVector(DN, 8, 1);

        data = [DeltLSBin; WNLSFBin; DNBin];
        word9 = gpsLNavWordEnc(data, wordNum);

        wordNum = 10;

        DeltLSF = utc.FutureLeapSecondCount;
        DeltLSFBin = HelperGPSNum2BitsVector(DeltLSF, 8, 1);

        data = [DeltLSFBin; zeros(14, 1)];
        word10 = gpsLNavWordEnc(data, wordNum);

        words = [word3; word4; word5; word6; word7; word8; word9; word10];
    case 25 % A-S flags/ SV configurations for 32 SVs, plus SV health for SV 25 through 32

        svid = subframe4SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);

        for iSV = 1:32 % Generate SVConfig for each PRN
            almdataidx = almanac.SatellitesPRNIndices==iSV;
            if nnz(almdataidx)
                svconfig = almanac.Data(almdataidx).SVConfig;
            else
                svconfig = 15;
            end
            svconfigBin = HelperGPSNum2BitsVector(svconfig, 4, 1); %#ok<NASGU>
            varName = ['svconfigBin' num2str(iSV)];
            eval([varName '=svconfigBin;']);
        end

        for iSV = 25:32 % Generate SVHealth for each PRN
            almdataidx = almanac.SatellitesPRNIndices==iSV;
            if nnz(almdataidx)
                svhealth = almanac.Data(almdataidx).SVHealth;
            else
                svhealth = 63;
            end
            svhealthBin = HelperGPSNum2BitsVector(svhealth, 6, 1); %#ok<NASGU>
            varName = ['svhealthBin' num2str(iSV)];
            eval([varName '=svhealthBin;']);
        end

        wordNum = 3;

        data = [dataid; svidBin; svconfigBin1; svconfigBin2; svconfigBin3; svconfigBin4];
        word3 = gpsLNavWordEnc(data, wordNum);

        for iword = 4:7
            varNum = (iword-4)*6+5:(iword-3)*6+4;
            eval(['data = [svconfigBin' num2str(varNum(1)) ';svconfigBin' num2str(varNum(2))...
                ';svconfigBin' num2str(varNum(3)) ';svconfigBin' num2str(varNum(4))...
                ';svconfigBin' num2str(varNum(5)) ';svconfigBin' num2str(varNum(6)) '];']);
            eval(['word' num2str(iword) '=gpsLNavWordEnc(data, iword);']);
        end

        wordNum = 8;

        data = [svconfigBin29; svconfigBin30; svconfigBin31; svconfigBin32;...
            zeros(2,1); svhealthBin25];
        word8 = gpsLNavWordEnc(data, wordNum);

        wordNum = 9;

        data = [svhealthBin26; svhealthBin27; svhealthBin28; svhealthBin29];
        word9 = gpsLNavWordEnc(data, wordNum);

        wordNum = 10;

        data = [svhealthBin30; svhealthBin31; svhealthBin32; zeros(4, 1)];
        word10 = gpsLNavWordEnc(data, wordNum);

        words = [word3; word4; word5; word6; word7; word8; word9; word10]; %#ok
end

subframe = [word1; word2; words];
end

function subframe = subframe5(opts)
% Generate 300 bits of GPS LNAV subframe5 data from the parameters given in
% 'opts' and 'almanac' input structures

almanac = opts.Almanac;

word1 = TLMWord(opts); % Having the default values of TLM word for speed
word2 = HandoverWord(opts);

dataid = [0; 1];
subframe5SVID = [(1:24)'; 51];

caseName = 'Empty';
if ismember(opts.PageID, 1:24)
    caseName = 'almanac';
elseif opts.PageID==25
    caseName = 'SVData';
end

switch caseName
    case 'almanac' % Almanac data for SV 25 through 32
        prnidx = subframe5SVID(opts.PageID);
        words = almstruct2bits(almanac, prnidx);
    case 'SVData'
        wordNum = 3;

        svid = subframe5SVID(opts.PageID);
        svidBin = HelperGPSNum2BitsVector(svid, 6, 1);

        toa = almanac.ReferenceTimeOfAlmanac;
        toaBin = HelperGPSNum2BitsVector(toa, 8, 2^12);

        WNa = mod(almanac.WeekNumModulo1024, 2^8);
        WNaBin = HelperGPSNum2BitsVector(WNa, 8, 1);

        data = [dataid; svidBin; toaBin; WNaBin];
        word3 = gpsLNavWordEnc(data, wordNum);

        for iSV = 1:24 % Generate SVHealth for each PRN
            almdataidx = almanac.SatellitesPRNIndices==iSV;
            if nnz(almdataidx)
                svhealth = almanac.Data(almdataidx).SVHealth;
            else
                svhealth = 63;
            end
            svhealthBin = HelperGPSNum2BitsVector(svhealth, 6, 1); %#ok<NASGU>
            varName = ['svhealthBin' num2str(iSV)];
            eval([varName '=svhealthBin;']);
        end

        for iword = 4:9
            varNum = (iword-4)*4+1:(iword-3)*4;
            eval(['data = [svhealthBin' num2str(varNum(1)) ';svhealthBin' num2str(varNum(2))...
                ';svhealthBin' num2str(varNum(3)) ';svhealthBin' num2str(varNum(4)) '];']);
            eval(['word' num2str(iword) '=gpsLNavWordEnc(data, iword);']);
        end

        wordNum =  10;
        data = zeros(22, 1);
        word10 = gpsLNavWordEnc(data, wordNum);

        words = [word3; word4; word5; word6; word7; word8; word9; word10];
end

subframe = [word1; word2; words];

end

function bits = almstruct2bits(almanac, PRNID)
%Converts one PRNID almanac data to bits from structure
%   BITS = almstruct2bits(ALMANAC, PRNID) converts almanac which is in a
%   structure to bits in word 3 through 10 along with proper parity for
%   LNAV data. This is useful in both subframe 4 and subframe 5.
%
% Parameters in ALMANAC
%   SatellitesPRNIndices
%   WeekNumModulo1024
%   ReferenceTimeOfAlmanac
%   Data.PRNID
%   Data.SVID
%   Data.URAID
%   Data.Eccentricity
%   Data.InclinationDifference
%   Data.RateOfRightAscension
%   Data.RootOfSemiMajorAxis
%   Data.LongitudeOfAscendingNode
%   Data.ArgumentOfPerigee
%   Data.MeanAnomaly
%   Data.SVClockCorrectionCoefficients
%   Data.SVHealth % 6 bits with MSB value to be health status flag
%   Data.NAVDataHealth % 3 bit NAV data health indications. If not specified,
%                      % default value of 3 three zero bits is taken
%   Data.SVConfig

wordNum = 3;

dataid = [0; 1];
almdataidx = almanac.SatellitesPRNIndices==PRNID;
if nnz(almdataidx)==0
    % Where PRNID is not found, then that PRNID satellite is a dummy
    % satellite and initialize the almanac message with alternating ones
    % and zeros and have the SV ID as zero.
    svidBin = zeros(6, 1);

    data = [dataid; svidBin; repmat([1; 0], 8, 1)];
    word3 = gpsLNavWordEnc(data,wordNum);

    wordNum = 4;
    data = repmat([1; 0], 12, 1);
    word4 = gpsLNavWordEnc(data,wordNum);

    wordNum = 5;
    word5 = gpsLNavWordEnc(data,wordNum);

    wordNum = 6;
    word6 = gpsLNavWordEnc(data,wordNum);

    wordNum = 7;
    word7 = gpsLNavWordEnc(data,wordNum);

    wordNum = 8;
    word8 = gpsLNavWordEnc(data,wordNum);

    wordNum = 9;
    word9 = gpsLNavWordEnc(data,wordNum);

    wordNum = 10;
    data = repmat([1; 0], 11, 1);
    word10 = gpsLNavWordEnc(data,wordNum);
else
    svidBin = HelperGPSNum2BitsVector(almanac.Data(almdataidx).PRNID, 6, 1);

    e = almanac.Data(almdataidx).Eccentricity;
    eBin = HelperGPSNum2BitsVector(e, 16, 2^-21);

    data = [dataid; svidBin; eBin];
    word3 = gpsLNavWordEnc(data,wordNum);

    wordNum = 4;

    toa = almanac.ReferenceTimeOfAlmanac;
    toaBin = HelperGPSNum2BitsVector(toa, 8, 2^(12));

    deli = almanac.Data(almdataidx).InclinationDifference; % delta_i is a correction to inclination term with reference to 0.3semicircles
    deliBin = HelperGPSNum2BitsVector(deli, 16, 2^(-19));

    data = [toaBin; deliBin];
    word4 = gpsLNavWordEnc(data,wordNum);

    wordNum = 5;

    OmegaDot = almanac.Data(almdataidx).RateOfRightAscension;
    OmegaDotBin = HelperGPSNum2BitsVector(OmegaDot, 16, 2^(-38));

    % SVHealth in almanac should be of 8 bits but SEM only gives
    % information of 5 bits. Other three bits is data health which is not
    % specified in SEM almanac format. So, if a user specifies data health
    % value, it will be taken. Else, data health will be three zero bits.
    if isfield(almanac.Data, 'NAVDataHealth')
        datahealth = almanac.Data(almdataidx).NAVDataHealth;
        datahealthBin = HelperGPSNum2BitsVector(datahealth, 3, 1);
    else
        datahealthBin = zeros(3, 1);
    end

    svhealth = almanac.Data(almdataidx).SVHealth;
    svhealthBin = HelperGPSNum2BitsVector(svhealth, 6, 1);

    data = [OmegaDotBin; datahealthBin; svhealthBin(2:end)];
    word5 = gpsLNavWordEnc(data,wordNum);

    wordNum = 6;

    sqrtA = almanac.Data(almdataidx).RootOfSemiMajorAxis;
    sqrtABin = HelperGPSNum2BitsVector(sqrtA, 24, 2^(-11));

    word6 = gpsLNavWordEnc(sqrtABin, wordNum);

    wordNum = 7;

    Omega_0 = almanac.Data(almdataidx).LongitudeOfAscendingNode;
    Omega_0Bin = HelperGPSNum2BitsVector(Omega_0, 24, 2^-23);

    word7 = gpsLNavWordEnc(Omega_0Bin, wordNum);

    wordNum = 8;

    omega = almanac.Data(almdataidx).ArgumentOfPerigee;
    omegaBin = HelperGPSNum2BitsVector(omega, 24, 2^-23);

    word8 = gpsLNavWordEnc(omegaBin, wordNum);

    wordNum = 9;

    M0 = almanac.Data(almdataidx).MeanAnomaly;
    M0Bin = HelperGPSNum2BitsVector(M0, 24, 2^-23);

    word9 = gpsLNavWordEnc(M0Bin, wordNum);

    wordNum = 10;

    af0 = almanac.Data(almdataidx).SVClockCorrectionCoefficients(1);
    af0Bin = HelperGPSNum2BitsVector(af0, 11, 2^-20);

    af1 = almanac.Data(almdataidx).SVClockCorrectionCoefficients(2);
    af1Bin = HelperGPSNum2BitsVector(af1, 11, 2^-38);

    data = [af0Bin(1:8); af1Bin; af0Bin(9:11)];
    word10 = gpsLNavWordEnc(data, wordNum);
end

bits = [word3; word4; word5; word6; word7; word8; word9; word10];
end

function word = TLMWord(tlmdata)
%HTLMWORD generates the telemetry (TLM) word
%   WORD = TLMWord(TLMDATA) generates the 30 bit length telemetry word,
%   WORD as per the data contained in TLMDATA. When TLMDATA.DefaultTLMFlag
%   is set, irrespective of anything in TLMWORD, default value of TLM word
%   will be generated. Leaving the 6 bits parity, default TLM word is
%   [1;0;0;0;1;0;1;1; zeros(16,1)]. TLMDATA is a structure or configuration
%   object with following parameters
%
%   TLMPreamble                    - an integer which is of 8 bit length
%   TLMMessage                     - an integer which is of 14 bits length
%   CEIDataSet.IntegrityStatusFlag - Binary value

preamble = HelperGPSNum2BitsVector(tlmdata.Preamble, 8, 1);
TLMMessage = HelperGPSNum2BitsVector(tlmdata.TLMMessage, 14, 1);
flags = [tlmdata.CEIDataSet.IntegrityStatusFlag; 0];
data = [preamble; TLMMessage; flags];
wordNumber = 1;
word = gpsLNavWordEnc(data,wordNumber);
end

function nmctBin = NMCT2Bin(erd)
%Converts the ERD values to binary values
%   NMCTBIN = NMCT2Bin(ERD) converts the ERD values to binary values and
%   returns a binary matrix NMCTBIN. ERD is a column vector of length 30
%   containing NMCT ERD values in each element of the array. NMCTBIN in a
%   binary matrix with number of rows same as that of ERD and number of
%   columns equal to 6.

numeleInNMCT = length(erd);

numBits = 6;
scaleFactor = 0.3;
nmctBin = zeros(numeleInNMCT, numBits);

for iERD = 1:numeleInNMCT
    nmctBin(iERD, :) = HelperGPSNum2BitsVector(erd(iERD), numBits, scaleFactor);
end
end

function codeWord = gpsLNavWordEnc(data,wordNumber)
%Encode the GPS LNAV word
%   CODEWORD = gpsLNavWordEnc(DATA,WORDNUMBER) encodes DATA using the
%   encoding scheme specified in IS-GPS-200K to obtain the encoded
%   CODEWORD. DATA is a vector of bits. If WORDNUMBER is either 2 or 10,
%   then DATA should be a vector of 22 bits and for other value of
%   WORDNUMBER, DATA is a vector of 24 bits. CODEWORD is a vector of 30
%   bits.

persistent d29Star d30Star
if wordNumber==1 || isempty(d29Star)
    d29Star=0;
    d30Star=0;
end
% Parity Polynomials for bits 25 to 30
d25 = [1 2 3 5 6 10 11 12 13 14 17 18 20 23];
d26 = [2 3 4 6 7 11 12 13 14 15 18 19 21 24];
d27 = [1 3 4 5 7 8 12 13 14 15 16 19 20 22];
d28 = [2 4 5 6 8 9 13 14 15 16 17 20 21 23];
d29 = [1 3 5 6 7 9 10 14 15 16 17 18 21 22 24];
d30 = [3 5 6 8 9 10 11 13 15 19 22 23 24];

data = data(:); % As data should always be a column vector

if wordNumber~=2 && wordNumber~=10
    %Vectorise parity application to data
    dataParity = d30Star*ones(24, 1);

    %Generate Codeword based on the parity equations from the ICD -
    %Hamming(32,26)
    codeWord = mod([dataParity+data; d29Star+sum(data(d25));...
        d30Star+sum(data(d26)); d29Star+sum(data(d27));...
        d30Star+sum(data(d28)); d30Star+sum(data(d29));...
        d29Star+sum(data(d30))],2);

    %Update the D30* and D29* Value for the next word of subframe
    d30Star = codeWord(30);
    d29Star = codeWord(29);

end

%Parity for the Handover word and word number 10 is generated differently
% Word number 2 and 10 have bits 29 and 30 set to 0 and the parity is then
% generated accordingly. The values of bits 23 and 24 are then calculated
% to satisfy this condition.
if wordNumber==2 || wordNumber==10
    %Vectorise Parity application to data
    dataParity = d30Star*ones(22, 1);

    %Bits 23 and 24 have to be calculated by setting D29 and D30 to zero
    bit24 = mod(sum(data(d29(d29<23)))+...
        d30Star,2);
    bit23 = mod(sum(data(d30(d30<23)))+bit24+...
        d29Star,2);

    %Generate Parity Encoded Data
    codeWord = mod([dataParity+data(1:22); d30Star+bit23;...
        d30Star+bit24; d29Star+sum(data(d25(d25<23)))+bit23;...
        d30Star+sum(data(d26(d26<23)))+bit24; d29Star+sum(data(d27(d27<23)));...
        d30Star+sum(data(d28(d28<23)))+bit23; 0; 0],2);

    %Update the D30* and D29* Value for the next word of subframe
    d29Star=0;
    d30Star=0;

end
end

function rem = modlocal(num, den)
% When NUM is exactly divisible by DEN, then output is DEN rather than 0.
% This is highly useful in applications where indexing cyclically rotates

rem = mod(num, den);
if rem == 0 && num ~= 0
    rem = den;
end
end

function CNAVBits = GPSCNAVDataEncode(cfg)
%GPSCNAVDataEncode Navigation data bits generation for L2C CNAV GPS
%
%   CNAVBITS = GPSCNAVDataEncode(CFG) generates the encoded bit
%   stream CNAVBITS which contains data bits for L2C GPS transmission.
%   CFG is a configuration object of type <a href="matlab:help('HelperGPSNavigationConfig')">HelperGPSNavigationConfig</a>. The
%   data bits are generated from the fields that go into modernized GPS
%   data. Each of these fields are available in the configuration object,
%   CFG. The properties and the encoding procedure is given in
%   Appendix III of IS-GPS-200L[1].

sequence = cfg.MessageTypes(:);
% Each message type contains 300 bits
numSequence = length(sequence);
CNAVBits = zeros(300*numSequence, 1, 'int8');

% Clearing persistent variables
getReducedAlmanacPacket();
getMidiAlmanacParams();
getCDCPacket();
getEDCPacket();

TOW = cfg.HOWTOW;
for idx = 1:numSequence
    % For every 6 seconds, 17 MSBs of TOW gets incremented by 1. As each
    % message is 12 seconds long, TOW must be incremented by 2.
    cfg.HOWTOW = TOW + (idx-1)*2;
    startIdx = 300*(idx - 1) + 1;
    endIdx = 300*(idx - 1) + 300;
    switch sequence(idx)
        case 10
            CNAVBits(startIdx:endIdx) = messageType10(cfg);
        case 11
            CNAVBits(startIdx:endIdx) = messageType11(cfg);
        case 30
            CNAVBits(startIdx:endIdx) = messageType30(cfg);
        case 31
            CNAVBits(startIdx:endIdx) = messageType31(cfg);
        case 32
            CNAVBits(startIdx:endIdx) = messageType32(cfg);
        case 33
            CNAVBits(startIdx:endIdx) = messageType33(cfg);
        case 34
            CNAVBits(startIdx:endIdx) = messageType34(cfg);
        case 35
            CNAVBits(startIdx:endIdx) = messageType35(cfg);
        case 36
            CNAVBits(startIdx:endIdx) = messageType36(cfg);
        case 37
            CNAVBits(startIdx:endIdx) = messageType37(cfg);
        case 12
            CNAVBits(startIdx:endIdx) = messageType12(cfg);
        case 13
            CNAVBits(startIdx:endIdx) = messageType13(cfg);
        case 14
            CNAVBits(startIdx:endIdx) = messageType14(cfg);
        case 15
            CNAVBits(startIdx:endIdx) = messageType15(cfg);
        otherwise
            CNAVBits(startIdx:endIdx) = messageType0(cfg);
    end
end
end

function bits = messageType0(opts)
% Generates bits for message type 0

bits = getCommonParamBits(opts, 0);
x = [ones(1,119);zeros(1,119)];
bits = [bits;x(:)];
bits = CRCGen(bits);
end

function bits = messageType10(opts)
% Generates bits for message type 10

bits = getCommonParamBits(opts, 10);
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.WeekNumber, 13, 1)];
bits = [bits;opts.CEIDataSet.SignalHealth(1)];
bits = [bits;opts.CEIDataSet.SignalHealth(2)];
bits = [bits;opts.CEIDataSet.SignalHealth(3)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceTimeCEIPropagation, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.URAEDID, 5, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceTimeOfEphemeris, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.SemiMajorAxisLength-26559710, 26, 2^-9)]; % What is transmitted is semi-major axis length difference with respect to ARef = 26559710
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ChangeRateInSemiMajorAxis, 25, 2^-21)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.MeanMotionDifference, 17, 2^-44)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.RateOfMeanMotionDifference, 23, 2^-57)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.MeanAnomaly, 33, 2^-32)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.Eccentricity, 33, 2^-34)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ArgumentOfPerigee, 33, 2^-32)];
bits = [bits;opts.CEIDataSet.IntegrityStatusFlag];
bits = [bits;opts.L2CPhasing];
bits = [bits;zeros(3,1)];
bits = CRCGen(bits);
end

function bits = messageType11(opts)
% Generates bits for message type 11

bits = getCommonParamBits(opts, 11);
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceTimeOfEphemeris, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.LongitudeOfAscendingNode, 33, 2^-32)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.Inclination, 33, 2^-32)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.RateOfRightAscension+2.6e-9, 17, 2^-44)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.InclinationRate, 15, 2^-44)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(1), 16, 2^-30)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(2), 16, 2^-30)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(3), 24, 2^-8)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(4), 24, 2^-8)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(5), 21, 2^-30)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.HarmonicCorrectionTerms(6), 21, 2^-30)];
bits = [bits;zeros(7,1)];
bits = CRCGen(bits);
end

function bits = messageType30(opts)
% Generates bits for message type 30

bits = getCommonParamBits(opts, 30);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.GroupDelayDifferential, 13, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.InterSignalCorrection(1), 13, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.InterSignalCorrection(2), 13, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.InterSignalCorrection(3), 13, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.InterSignalCorrection(4), 13, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Alpha(1), 8, 2^-30)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Alpha(2), 8, 2^-27)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Alpha(3), 8, 2^-24)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Alpha(4), 8, 2^-24)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Beta(1), 8, 2^11)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Beta(2), 8, 2^14)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Beta(3), 8, 2^16)];
bits = [bits;HelperGPSNum2BitsVector(opts.Ionosphere.Beta(4), 8, 2^16)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceWeekNumberCEIPropagation, 8, 1)];
bits = [bits;zeros(12,1)];
bits = CRCGen(bits);
end

function bits = messageType31(opts)
% Generates bits for message type 31

Almanac = opts.Almanac;
bits = getCommonParamBits(opts, 31);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(mod(Almanac.WeekNumModulo1024, 2^14), 13, 1)];
bits = [bits;HelperGPSNum2BitsVector(Almanac.ReferenceTimeOfAlmanac, 8, 2^12)];
prevDummySat = 0;
reducedAlmanacPackets = zeros(124, 1);
for idx = 1:4
    startIdx = (idx-1) * 31 + 1;
    endIdx = (idx - 1) * 31 + 31;
    [reducedAlmanacBits, dummySat] = getReducedAlmanacPacket(Almanac, prevDummySat);
    reducedAlmanacPackets(startIdx:endIdx) = reducedAlmanacBits;
    prevDummySat = dummySat;
end
bits = [bits;reducedAlmanacPackets];
bits = [bits;zeros(4,1)];
bits = CRCGen(bits);
end

function bits = messageType32(opts)
% Generates bits for message type 32

bits = getCommonParamBits(opts, 32);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.ReferenceTimeEOP, 16, 2^4)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.XAxisPolarMotionValue, 21, 2^-20)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.XAxisPolarMotionDrift, 15, 2^-21)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.YAxisPolarMotionValue, 21, 2^-20)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.YAxisPolarMotionDrift, 15, 2^-21)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.UT1_UTCDifference, 31, 2^-23)];
bits = [bits;HelperGPSNum2BitsVector(opts.EarthOrientation.RateOfUT1_UTCDifference, 19, 2^-25)];
bits = [bits;zeros(11,1)];
bits = CRCGen(bits);
end


function bits = messageType33(opts)
% Generates bits for message type 33

bits = getCommonParamBits(opts, 33);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.UTCTimeCoefficients(1), 16, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.UTCTimeCoefficients(2), 13, 2^-51)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.UTCTimeCoefficients(3), 7, 2^-68)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.PastLeapSecondCount, 8, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.ReferenceTimeUTCData, 16, 2^4)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.TimeDataReferenceWeekNumber, 13, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.LeapSecondReferenceWeekNumber, 13, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.LeapSecondReferenceDayNumber, 4, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.UTC.FutureLeapSecondCount, 8, 1)];
bits = [bits;zeros(51,1)];
bits = CRCGen(bits);
end

function bits = messageType34(opts)
% Generates bits for message type 34

DCParams = opts.DifferentialCorrection;
bits = getCommonParamBits(opts, 34);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCDataPredict, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCData, 11, 300)];
bits = [bits;getCDCPacket(DCParams)];
EDCPacket = getEDCPacket(DCParams);
bits = [bits;EDCPacket(2:end)];
bits = CRCGen(bits);
end

function bits = messageType35(opts)
% Generates bits for message type 35

bits = getCommonParamBits(opts, 35);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.ReferenceTimeGGTO, 16, 2^4)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.WeekNumberGGTO, 13, 2^0)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.GNSSID, 3, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.GGTOCoefficients(1), 16, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.GGTOCoefficients(2), 13, 2^-51)];
bits = [bits;HelperGPSNum2BitsVector(opts.TimeOffset.GGTOCoefficients(3), 7, 2^-68)];
bits = [bits;zeros(81,1)];
bits = CRCGen(bits);
end

function bits = messageType36(opts)
% Generates bits for message type 36

bits = getCommonParamBits(opts, 36);
bits = [bits;getClockParams(opts)];
numCharInMessage = 18;
fullText = char(opts.TextInMessageType36(:).');
if length(fullText) >= numCharInMessage
    text = fullText(1:numCharInMessage);
else
    text = [fullText, repmat(' ',1, numCharInMessage-length(fullText))];
end
bits = [bits;GPSText2Bits(text)];

% Text message and text page are two separate fields. IS-GPS-200k does not
% explain what this text page field specifies. Filled it as zeros of length
% 4 where 4 is the length of the field.
bits = [bits;zeros(4,1)];
bits = [bits;zeros(1,1)];
bits = CRCGen(bits);
end

function bits = messageType37(opts)
% Generates bits for message type 37

Almanac = opts.Almanac;
bits = getCommonParamBits(opts, 37);
bits = [bits;getClockParams(opts)];
bits = [bits;HelperGPSNum2BitsVector(mod(Almanac.WeekNumModulo1024, 2^14), 13, 1)];
bits = [bits;HelperGPSNum2BitsVector(Almanac.ReferenceTimeOfAlmanac, 8, 2^12)];
bits = [bits;getMidiAlmanacParams(Almanac)];
bits = CRCGen(bits);
end

function bits = messageType12(opts)
% Generates bits for message type 12

Almanac = opts.Almanac;
bits = getCommonParamBits(opts, 12);
bits = [bits;HelperGPSNum2BitsVector(mod(Almanac.WeekNumModulo1024, 2^14), 13, 1)];
bits = [bits;HelperGPSNum2BitsVector(Almanac.ReferenceTimeOfAlmanac, 8, 2^12)];
prevDummySat = 0;
reducedAlmanacPackets = zeros(217, 1);
for idx = 1:7
    startIdx = (idx-1) * 31 + 1;
    endIdx = (idx - 1) * 31 + 31;
    [reducedAlmanacBits, dummySat] = getReducedAlmanacPacket(Almanac, prevDummySat);
    reducedAlmanacPackets(startIdx:endIdx) = reducedAlmanacBits;
    prevDummySat = dummySat;
end
bits = [bits;reducedAlmanacPackets];
bits = CRCGen(bits);

end

function bits = messageType13(opts)
% Generates bits for message type 13

DCParams = opts.DifferentialCorrection;
bits = getCommonParamBits(opts, 13);
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCDataPredict, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCData, 11, 300)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;getCDCPacket(DCParams)];
bits = [bits;zeros(6,1)];
bits = CRCGen(bits);
end

function bits = messageType14(opts)
% Generates bits for message type 14

DCParams = opts.DifferentialCorrection;
bits = getCommonParamBits(opts, 14);
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCDataPredict, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(DCParams.ReferenceTimeDCData, 11, 300)];
bits = [bits;getEDCPacket(DCParams)];
bits = [bits;getEDCPacket(DCParams)];
bits = [bits;zeros(30,1)];
bits = CRCGen(bits);
end

function bits = messageType15(opts)
% Generates bits for message type 15

numCharInMessage = 29;
bits = getCommonParamBits(opts, 15);
fullText = char(opts.TextInMessageType15(:).');
if length(fullText) >= numCharInMessage
    text = fullText(1:numCharInMessage);
else
    text = [fullText, repmat(' ',1, numCharInMessage-length(fullText))];
end
bits = [bits;GPSText2Bits(text)];

% Text message and text page are two separate fields. IS-GPS-200k does not
% explain what this text page field specifies. Filled it as zeros of length
% 4 where 4 is the length of the field.
bits = [bits;zeros(4,1)];
bits = [bits;zeros(2,1)];
bits = CRCGen(bits);
end

function bits = getCommonParamBits(opts, messageTypeId)
% Generates bits corresponding to preamble, PRNIdx, message type ID, TOW
% and alert flag which are present in all message types

bits = HelperGPSNum2BitsVector(opts.Preamble, 8, 1);
bits = [bits;HelperGPSNum2BitsVector(opts.PRNID, 6, 1)];
bits = [bits;HelperGPSNum2BitsVector(messageTypeId, 6, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.HOWTOW, 17, 1)];
bits = [bits;opts.CEIDataSet.AlertFlag];
end

function bits = getClockParams(opts)
% Generates clock correction and accuracy parameters bits for message types
% 30 to 37.

bits = HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceTimeCEIPropagation, 11, 300);
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.URANEDID(1), 5, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.URANEDID(2), 3, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.URANEDID(3), 3, 1)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.ReferenceTimeOfClock, 11, 300)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.SVClockCorrectionCoefficients(1), 26, 2^-35)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.SVClockCorrectionCoefficients(2), 20, 2^-48)];
bits = [bits;HelperGPSNum2BitsVector(opts.CEIDataSet.SVClockCorrectionCoefficients(3), 10, 2^-60)];
end

function bits = GPSText2Bits(text)
TempBits = comm.internal.utilities.de2biBase2LeftMSB((double(text(:))),8)';
bits = TempBits(:);
end

function bits = CRCGen(inputBits)
% Returns 300-bit message after appending 24 CRC bits at the end of message

persistent crcgenerator
if isempty(crcgenerator)
    crcgenerator = comm.CRCGenerator('Polynomial', ...
        'z^24 + z^23 + z^18 + z^17 + z^14 + z^11 + z^10 + z^7 + z^6 + z^5 + z^4 + z^3 + z + 1');
end
bits = crcgenerator(inputBits);
end

function [bits, dummySat] = getReducedAlmanacPacket(Almanac, prevDummySat)
%GETREDUCEDALMANACPACKET generated bits in reduced almanac packet for
%   message types 31 and 12.
%   [BITS, DUMMYSAT] = getReducedAlmanac(ALMANAC, PREVDUMMYSAT) generates
%   31 bits reduced almanac packet for an SV.
%
%   BITS            - Reduced almanac packet
%   DUMMYSAT        - Flag indicating current packet contains
%                     dummy satellite data
%   ALMANAC         - Structure containing almanac parameters
%   PREVDUMMYSAT    - Flag indicating dummy satellite data in any of the
%                     previous reduced almanac packets in the message

% Index to transmit reduced almanac data in ascending order of PRN numbers
% of satellites in constellation
persistent redAlmanacIdx nextPRN lastDummyBit

if nargin == 0
    clear redAlmanacIdx nextPRN lastDummyBit
else
    if isempty(redAlmanacIdx) || redAlmanacIdx >= 31
        redAlmanacIdx = 0;
        nextPRN = 1;
    end
    if isempty(lastDummyBit)
        lastDummyBit = 1;
    end

    redAlmanacIdx = redAlmanacIdx + 1;
    PRN_a = Almanac.Data(redAlmanacIdx).PRNID;

    if prevDummySat == 1
        % If there is a previous dummy satellite in the message, bits after PRN
        % field of the reduced almanac packet of dummy satellite through the last
        % bit of last almanac packet will be filler bits.
        x = [xor(ones(1,15),lastDummyBit);repmat(lastDummyBit,1,15)];
        bits = [x(:);xor(1,lastDummyBit)];
        lastDummyBit = xor(1,lastDummyBit);
        dummySat = 1;
        redAlmanacIdx = redAlmanacIdx - 1;

    elseif PRN_a == nextPRN && prevDummySat == 0
        bits = HelperGPSNum2BitsVector(PRN_a, 6, 1);
        del_A = (Almanac.Data(redAlmanacIdx).RootOfSemiMajorAxis^2) - 26559710; % del_A is with respect to 26559710 meters
        bits = [bits;HelperGPSNum2BitsVector(del_A, 8, 2^9)];
        omega_0 = Almanac.Data(redAlmanacIdx).LongitudeOfAscendingNode;
        bits = [bits;HelperGPSNum2BitsVector(omega_0, 7, 2^-6)];
        phi_0 = Almanac.Data(redAlmanacIdx).MeanAnomaly + Almanac.Data(redAlmanacIdx).ArgumentOfPerigee;
        bits = [bits;HelperGPSNum2BitsVector(phi_0, 7, 2^-6)];
        bits = [bits;Almanac.Data(redAlmanacIdx).L1Health];
        bits = [bits;Almanac.Data(redAlmanacIdx).L2Health];
        bits = [bits;Almanac.Data(redAlmanacIdx).L5Health];
        dummySat = 0;
        nextPRN = nextPRN + 1;

    elseif PRN_a ~= nextPRN && prevDummySat == 0
        % If current satellite is a dummy satellite, PRN field will be "000000"
        % and all subsequent bits through the last bit of the last packet in
        % the message will be filler bits.
        bits = zeros(6, 1);
        x = [ones(1,12);zeros(1,12)];
        bits = [bits;x(:);1];
        dummySat = 1;
        lastDummyBit = 1;
        nextPRN = nextPRN + 1;
        redAlmanacIdx = redAlmanacIdx - 1;
    end
end
end

function bits = getMidiAlmanacParams(Almanac)
% Generates midi almanac parameters bits for message type 37.

% Index to transmit midi almanac data in ascending order of PRN numbers of
% satellites in constellation
persistent midiAlmanacIdx
if nargin == 0
    clear midiAlmanacIdx
else
    if isempty(midiAlmanacIdx) || midiAlmanacIdx >= 31
        midiAlmanacIdx = 0;
    end
    midiAlmanacIdx = midiAlmanacIdx + 1;

    PRN_a = Almanac.Data(midiAlmanacIdx).PRNID;
    bits = HelperGPSNum2BitsVector(PRN_a, 6, 1);
    bits = [bits;Almanac.Data(midiAlmanacIdx).L1Health];
    bits = [bits;Almanac.Data(midiAlmanacIdx).L2Health];
    bits = [bits;Almanac.Data(midiAlmanacIdx).L5Health];
    e = Almanac.Data(midiAlmanacIdx).Eccentricity;
    bits = [bits;HelperGPSNum2BitsVector(e, 11, 2^-16)];
    delta_i = Almanac.Data(midiAlmanacIdx).InclinationDifference;
    bits = [bits;HelperGPSNum2BitsVector(delta_i, 11, 2^-14)];
    omegaDot = Almanac.Data(midiAlmanacIdx).RateOfRightAscension;
    bits = [bits;HelperGPSNum2BitsVector(omegaDot, 11, 2^-33)];
    sqrtA = Almanac.Data(midiAlmanacIdx).RootOfSemiMajorAxis;
    bits = [bits;HelperGPSNum2BitsVector(sqrtA, 17, 2^-4)];
    omega_0 = Almanac.Data(midiAlmanacIdx).LongitudeOfAscendingNode;
    bits = [bits;HelperGPSNum2BitsVector(omega_0, 16, 2^-15)];
    omega = Almanac.Data(midiAlmanacIdx).ArgumentOfPerigee;
    bits = [bits;HelperGPSNum2BitsVector(omega, 16, 2^-15)];
    M_0 = Almanac.Data(midiAlmanacIdx).MeanAnomaly;
    bits = [bits;HelperGPSNum2BitsVector(M_0, 16, 2^-15)];
    a_f0 = Almanac.Data(midiAlmanacIdx).SVClockCorrectionCoefficients(1);
    bits = [bits;HelperGPSNum2BitsVector(a_f0, 11, 2^-20)];
    a_f1 = Almanac.Data(midiAlmanacIdx).SVClockCorrectionCoefficients(2);
    bits = [bits;HelperGPSNum2BitsVector(a_f1, 10, 2^-37)];
end
end

function bits = getCDCPacket(DCParams)
% Generate clock differential correction(CDC) packet for an SV

% Index to transmit CDC data in ascending order of PRN numbers of
% satellites in constellation
persistent CDCIdx

if nargin == 0
    clear CDCIdx
else
    if isempty(CDCIdx) || CDCIdx >= 31
        CDCIdx = 0;
    end

    CDCIdx = CDCIdx + 1;
    PRNIdx = DCParams.Data(CDCIdx).PRNID;
    if ~isempty(PRNIdx)
        bits = DCParams.Data(CDCIdx).DCDataType;
        bits = [bits;HelperGPSNum2BitsVector(PRNIdx, 8, 1)];
        delta_a_f0 = DCParams.Data(CDCIdx).SVClockBiasCoefficient;
        bits = [bits;HelperGPSNum2BitsVector(delta_a_f0, 13, 2^-35)];
        delta_a_f1 = DCParams.Data(CDCIdx).SVClockDriftCorrection;
        bits = [bits;HelperGPSNum2BitsVector(delta_a_f1, 8, 2^-51)];
        UDRA = DCParams.Data(CDCIdx).UDRAID;
        bits = [bits;HelperGPSNum2BitsVector(UDRA, 5, 1)];

    else
        % In case of dummy satellite, PRN field is filled with all 1s and the
        % subsequent bits till the end of packet are filler bits.
        bits = [0;ones(8, 1)];
        x = [ones(1,13);zeros(1,13)];
        bits = [bits;x(:)];
    end
end
end

function bits = getEDCPacket(DCParams)
% Generates Ephemeris differential correction(EDC) packet for an SV

% Index to transmit EDC data in ascending order of PRN numbers of
% satellites in constellation
persistent EDCIdx
if nargin == 0
    clear EDCIdx
else
    if isempty(EDCIdx) || EDCIdx >= 31
        EDCIdx = 0;
    end

    EDCIdx = EDCIdx + 1;
    PRNIdx = DCParams.Data(EDCIdx).PRNID;
    if ~isempty(PRNIdx)
        bits = DCParams.Data(EDCIdx).DCDataType;
        bits = [bits;HelperGPSNum2BitsVector(PRNIdx, 8, 1)];
        alphaCorrection = DCParams.Data(EDCIdx).AlphaCorrection;
        bits = [bits;HelperGPSNum2BitsVector(alphaCorrection, 14, 2^-34)];
        betaCorrection = DCParams.Data(EDCIdx).BetaCorrection;
        bits = [bits;HelperGPSNum2BitsVector(betaCorrection, 14, 2^-34)];
        gammaCorrection = DCParams.Data(EDCIdx).GammaCorrection;
        bits = [bits;HelperGPSNum2BitsVector(gammaCorrection, 15, 2^-32)];
        delta_i = DCParams.Data(EDCIdx).InclinationCorrection;
        bits = [bits;HelperGPSNum2BitsVector(delta_i, 12, 2^-32)];
        delta_omega = DCParams.Data(EDCIdx).RightAscensionCorrection;
        bits = [bits;HelperGPSNum2BitsVector(delta_omega, 12, 2^-32)];
        delta_A = DCParams.Data(EDCIdx).SemiMajorAxisCorrection;
        bits = [bits;HelperGPSNum2BitsVector(delta_A, 12, 2^-9)];
        UDRADot = DCParams.Data(EDCIdx).UDRARateID;
        bits = [bits;HelperGPSNum2BitsVector(UDRADot, 5, 1)];

    else
        % In case of dummy satellite, PRN field is filled with all 1s and the
        % subsequent bits till the end of packet are filler bits.
        bits = [0;ones(8, 1)];
        x = [ones(1,42);zeros(1,42)];
        bits = [bits;x(:)];
    end
end
end
